Imagery data is widespread and includes for example 3D seismic data, HS image, fingerprints and multimedia images. In particular, 3D seismic data and HS images are considered to be very large datasets. Efficient algorithms are needed to compress such datasets in order to store or transmit the data via different networks (wireless, Internet, TCP/IP, etc), or from various moving entities (e.g. planes, satellites, boats or cars) to other entities such as base stations. The compression ratio between the size of an original input image and a compressed image should be as high as possible without damaging the interpretation procedures applied after decompression. The compression of seismic, HS datasets and fingerprints should preserve fine details which are critical for interpretations after decompression. Interpretation is less critical for multimedia type images, where compression/decompression should however preserve acceptable visual quality. Multimedia images usually do not go through interpretation analysis after decompression.
A typical known data compression scheme is shown in FIG. 1a. This scheme operates on a 2D imagery data input. 2D data sets may be defined by two, first and second orthogonal axes or “directions”. Commonly, these axes are referred to as “horizontal” and “vertical”. The scheme includes two phases: 1) application of a transform to the horizontal direction of the image, step 100a, followed by the application of a transform to the vertical direction of the image, step 110a, both steps being part of a lossless phase. In general, the same transform (e.g. a wavelet transform) is applied to both directions; 2) processing of horizontal and vertical wavelet transform coefficients to obtain a 2D compressed image, step 120a. FIG. 1b gives more details of processing step 120a. Step 120a includes two phases: 1) quantization of the transform coefficients in a lossy phase step 130b; and 2) application of entropy coding for packing the quantized coefficients into a compressed form in a lossless phase step 140b. The output of step 140b is a 2D compressed image. Decompression is done in a reverse order.
Wavelet transforms have a successful record in achieving high compression ratios for still images while achieving high quality. For example, the transform used to perform steps 100a and 110a in FIG. 1a above is typically the 2D wavelet transform of the JPEG2000 standard which replaces the regular JPEG. The wavelet transform is the transform of choice in many compression algorithms. For example, the wavelet transform was chosen to handle seismic compression, see e.g. P. L. Donoho, R. A. Ergas, and J. D. Villasenor, “High-performance seismic trace compression”, in 65th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, (1995), 160-163; M. F. Khene and S. H. Abdul-Jauwad, “Efficient seismic compression using the lifting scheme”, in 70th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, (2000), 2052-2054; and A. Vassiliou and M. V. Wickerhauser, “Comparison of wavelet image coding schemes for seismic data compression”, in 67th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, (1997), 1334-1337.
The coding schemes of wavelet transforms described in J. M. Shapiro, “Embedded image coding using zerotree of wavelet coefficients”, IEEE Trans. Sign. Proc., vol. 41, pp. 3445-3462, 1993 (hereinafter “EZW”) and A. Said and W. W. Pearlman, “A new, fast and efficient image codec based on set partitioning in hierarchical trees” IEEE Trans. on Circ. and Syst. for Video Tech., vol. 6, pp. 243-250, 1996 (hereinafter “SPIHT”), are used in well known 2D wavelet transform schemes for still image compression. The wavelet transform in used in both directions. Both EZW and SPIHT utilize the correlation among the multiscale (also called multiresolution) decomposition of an image to achieve high compression ratios. The SPIHT scheme is associated with wavelet transforms that rely on the space-frequency localization of the wavelets and the tree structure of the coefficient arrays in its multiscale representation.
Wavelet transforms are also used to compress HS data cubes. A HS data cube can be captured by a special HS camera. The HS camera captures the same image (also called “spatial image”) in many wavebands. The number of wavebands depends on the camera resolution and can vary from a few wavebands (“multi-spectral imaging”) to a few thousands of wavebands (“HS imaging”). A HS image can be viewed as a 3D image cube where the X-Y plane is the spatial description of the image and the Z direction are the wavebands. In waveband scanning, which captures the HS spectra from the X-Z plane, a HS camera collects each time unit a line of all the wavebands (the Z direction). Each spatial pixel is represented by a vector (also called hereinafter multipixel) of the intensities in all the available wavebands. 2D compression of planes (wavebands×multipixels) is preferable since then no buffering is needed, because compression is done according to the data capturing mechanism. Compression of HS cubes should retain the spectral characteristic features of the multipixels. A number of HS compression algorithms were suggested such as G. Motta, F. Rizzo, and J. A. Storer (editors), “Hyperspectral data compression”, Kluwer Academic Publishers, 2006, 273-308. Many of them extend wavelet-based compression methods to compress 3D HS data cubes such as Y. H Tseng, H. K. Shih, and P. H. Hsu, “Hyperspectral Image Compression Using Three-dimensional Wavelet Transformation”, Proceedings of the 21th Asia Conference on Remote Sensing, 2000, pp. 809-814, X. Tang and W. A. Pearlman, “Three-dimensional wavelet-based compression of hyperspectral images”, in “Hyperspectral data compression”, Eds. G. Motta, F. Rizzo, and J. A. Storer, Kluwer Acad. Publ., 2006, 273-308, J. E. Fowler, J. T. Rucker, “Three-Dimensional Wavelet-Based Compression of Hyperspectral Imagery”, in “Hyperspectral Data Exploitation”, Ed. Chein-I Chang, John Wiley & Sons, (2007), 379-407.
Due to the great variability of seismic, HS and fingerprints images, which have inherently noisy backgrounds, oscillatory nature and many fine details that have to be preserved, wavelet based methods do not produce high compression ratios, see e.g. A. Z. Averbuch, F. Meyer, J-O. Stromberg, R. Coifman and A. Vassiliou, “Efficient Compression for Seismic Data”, IEEE Trans. on Image Processing, vol. 10, no. 12, pp. 1801-1814, 2001 (hereinafter “AMSCV”) and F. G. Meyer, “Fast compression of seismic data with local trigonometric bases”, Proc. SPIE 3813, Wavelet Applications in Signal and Image Processing VII, A. Aldroubi, A. F. Laine, and M. A. Unser, Eds., 1999, pp. 648-658 (hereinafter “FM”). Both AMSCV and FM showed that wavelet transforms do not fit seismic compression because of the oscillatory patterns present in seismic data, and that wavelet transforms do not capture efficiently these patterns. During a compression process, these patterns necessitate many bits to preserve features which are critical for their interpretation, thus no high compression ratios are achieved.
Seismic data has different structure in its two directions. While in one direction (say the horizontal one) the structure is piece-wise smooth and thus suits the capability of a wavelet transform to make sparse piece-wise smooth data, the traces in the other (vertical) direction include oscillatory patterns. While wavelets provide a sparse representation in the horizontal direction, helping to achieve a high compression ratio of the seismic data, they fail to properly (properly means here fewer coefficients) represent the vertical oscillations. The same is true for HS data cubes and fingerprint images.
On the other hand, the local cosine transform or “LCT” (R. R. Coifman and Y. Meyer, “Remarques sur l'analyse de Fourier a fenetre”, C. R. Acad. Sci., pp. 259-261, 1991 (hereinafter “CMLCT”) catches well oscillatory patterns. 2D LCT (in both horizontal and vertical directions) has been applied for seismic compression, see e.g. AMSCV, FM, Y. Wang and R.-S. Wu, “Seismic data compression by an adaptive local cosine/sine transform and its effect on migration”, Geophysical Prospecting, vol. 48, pp. 1009-1031, 2000 (hereinafter WWU) and V. A. Zheludev, D. D. Kosloff, E. Y. Ragoza, “Compression of segmented 3D seismic data”, International Journal of Wavelets, Multiresolution and Information Processing, vol. 2, no. 3, (2004), 269-281. The following two discrete cosine transforms (hereinafter “DCT”) types are used in image compression (see K. R. Rao and P. Yip, Discrete Cosine Transform, Academic press, New York, 1990):                1. DCT-II of the signal f={fn}n=0N−1, N=2p, used as the transform step (steps 100a and 110a in FIG. 1a) in the JPEG standard for still image compression, as        
                    f        II            ⁡              (        k        )              =                  b        ⁡                  (          k          )                    ⁢                        ∑                      n            =            0                                N            -            1                          ⁢                              f            n                    ⁢                      cos            ⁡                          [                                                                    k                    ⁢                                                                                  ⁢                    π                                    N                                ⁢                                  (                                      n                    +                                          1                      2                                                        )                                            ]                                            ,          ⁢            b      ⁡              (        k        )              =          {                                                                                    1                  /                                      2                                                                                                                    if                    ⁢                                                                                  ⁢                    k                                    =                  0                                                                                    1                                            otherwise                                              ⁢                                          ⁢                      f            n            II                          =                              2            N                    ⁢                                    ∑                              k                =                0                                            N                -                1                                      ⁢                                          b                ⁡                                  (                  k                  )                                            ⁢                                                f                  II                                ⁡                                  (                  k                  )                                            ⁢                                                cos                  ⁡                                      [                                                                                            k                          ⁢                                                                                                          ⁢                          π                                                N                                            ⁢                                              (                                                  n                          +                                                      1                            2                                                                          )                                                              ]                                                  .                                                                        2. DCT-IV of the signal f={fn}n=0N−1, N=2p, used as the transform step (steps 100a and 110a in FIG. 1a) in LCT as        
                    f        IV            ⁡              (        k        )              =                  ∑                  n          =          0                          N          -          1                    ⁢                        f          n                ⁢                  cos          ⁡                      [                                          π                N                            ⁢                              (                                  k                  +                                      1                    2                                                  )                            ⁢                              (                                  n                  +                                      1                    2                                                  )                                      ]                                ,          ⁢            f      n      IV        =                  2        N            ⁢                        ∑                      k            =            0                                N            -            1                          ⁢                                            f              IV                        ⁡                          (              k              )                                ⁢                                    cos              ⁡                              [                                                      π                    N                                    ⁢                                      (                                          k                      +                                              1                        2                                                              )                                    ⁢                                      (                                          n                      +                                              1                        2                                                              )                                                  ]                                      .                              DCT-IV is a good choice for coding oscillatory signals. The basis functions of DCT-IV are even on the left side with respect to
  -      1    2  and odd on the right side with respect to
  N  -            1      2        .  Therefore, direct application of the DCT-IV to partitioned data leads to severe boundary discrepancies. However, this transform serves as a base for LCTs which are window lapped DCT-IV DCTs. These bases were successfully exploited for image compression in general and seismic data in particular, see e.g. AMSCV, FM, WWU, A. Averbuch, G. Aharoni, R. Coifman and M. Israeli, “Local cosine transform—A method for the reduction of blocking effects in JPEG”, Journal of Mathematical Imaging and Vision, Special Issue on Wavelets, Vol. 3, pp. 7-38, 1993 (hereinafter AJPEG) and G. Matviyenko, “Optimized local trigonometric bases”, Applied and Computational Harmonic Analysis, 3, 301-323, 1996, (hereinafter MAT).
Assume we have a signal s={sk}k=0N−1 and some partition P of the interval 0: N−1. The idea behind the lapped DCT-IV is to apply overlapped bells to adjacent sub-intervals. Then, the overlapping parts are folded back to the sub-intervals across the endpoints of the sub-intervals and the DCT-IV on each sub-interval is implemented. In the reconstruction phase, the transform coefficients are unfolded. For details, see CMLCT and AMSCV. This transform is called P-based LCT. There are many available bells. Their descriptions are given in CMLCT, AJPEG and MAT. As an example, we use in our experiments the bell
      b    ⁡          (      x      )        =      sin    ⁢          π      2        ⁢                  (                  x          +                      1            2                          )            .      Wavelet Transforms
The multiscale wavelet transform of a signal is implemented via iterated multirate filtering. One step in the transform of a signal S={sk}k=0N−1 of length N consists of filtering the signal using a half-band low-pass filter L and a half-band high-pass filter H, followed by factor 2 down-sampling of both filtered signals. This produces two blocks of coefficients wL1={lk1}k=0N/2−1 and wH1={hk1}k=0N/2−1, each of length N/2. Coefficients wL1 include the entire information of the low frequency component of signal S while coefficients wH1 do the same for the high frequency component. Blocks wL1 and wH1 cut the Nyquist frequency band F of the signal S:F→FL1∪FH1 into half.
The wavelet transform coefficients also have a spatial meaning. The coefficient lm1 is the result of weighted averaging of the set Sm1,λ of signal samples. The set Sm1,λ is centered around sample s2m of signal S and its scope λ1 is equal to the width of the impulse response (IR) of filter L, provided L is a finite impulse response (FIR) filter. If L is an infinite impulse response (IIR) filter whose IR decays rapidly, then λ equals the effective width of the IR of filter L. The coefficient hm1 is the result of numerical differentiation of some order d of signal S at point 2m+1. For this, the set Sm1,χ of signal samples is involved, whose scope χ1 is equal to the (effective) width of the IR of filter H.
The next step of the wavelet transform applies the pair of filters L and H to the coefficients array wL1. The filtering is followed by down-sampling. The produced blocks of coefficients wL2={lk2} and wH2={hk2} from, respectively, L and H, are of length N/4 each, and cut in half the sub-band FL1→FL2∪FH2FL→FL2∪FH2∪FH1. In the time domain, coefficient lm2 is the result from averaging the set of samples Sm2,λ. The set Sm2,λ is centered around sample s4m of signal S and its scope is λ2≈2λ1. Similarly, coefficient hm2 is associated with the set Sm2,χ of samples, which is centered around sample s4m+2 of signal S and its scope is χ2≈2χ1. Note that set Sm2,λ occupies approximately the same area as the pair of its “offspring” sets S2m1,χ and S2m+11,χ.
Next, this decomposition is iterated to reach scale J. It produces the coefficients array wJ whose structure iswJ=wLJ∪wHJ∪wHJ−1∪ . . . ∪wH1.  (1)Respectively, the Nyquist frequency band F is split into sub-bands whose widths are distributed in a logarithmic way to becomeF→FLJ∪FHJ∪FHJ−1∪ . . . ∪FH1.  (2)The diagram of a three-scale wavelet transform and the layout of the transform coefficients are displayed in FIG. 2.
The wavelet transform of a 2D array T={tn,m} of size N×M is implemented in a tensor product. In other words, the 1D wavelet transform is applied to each side (direction). First, the pair of filters L and H is applied to the columns of T and the results are down-sampled. This yields coefficients arrays L and H of size N/2×M. Then, filters L and H are applied to the rows of L and H. This filtering is followed by down-sampling which results in four sub-array coefficients LL, LH, HL, HH of size N/2×M/2. The 2D Nyquist frequency domain is split accordingly. Then, the above procedure is applied to the coefficient array LL to produce the sub-arrays (LL)LL, (LL)LH, (LL)HL, (LL)HH of size N/4×M/4. Then, this procedure is iterated using (LL)LL instead of LL and so on. The layout of the transform coefficients, which corresponds to the Nyquist frequency partition, for the three-scale wavelet transform, is displayed in FIG. 3.
SPIHT Coding
The coefficients from (LL)LH, LL(HL), LL(HH) as shown in FIG. 3 are produced by the same operations as, respectively, the coefficients of LH, HL, HH. The coefficient cnmllhlε(LL)HL from the second scale is associated with the set Pnmllhl of pixels which is centered around pixel p4n+2,4m. Set Pnmllhl occupies (approximately) the same area as the four sets Pυ,μhl which are associated with the first scale coefficients cυ,μhlεHL, υ=2n+1,2n+3, μ=2m−2,2m+2, occupy. In this sense, cυ,μhlεHL are the descendants of coefficient cnmllhlε(LL)HL. Similar relations exist between the coefficients of LH, HH and (LL)LH, (LL)HH and also between the coefficients of other adjacent scales. These relations are illustrated in FIG. 4.
The ancestor-descendant relationship between wavelet transform coefficients in different scales, where the coefficients are located at the same spatial area as shown in FIG. 4, is exploited in the EZW codec. The EZW codec takes advantage of the self-similarity between wavelet coefficients across the decomposed scales and their decay toward high frequency. The EZW codec relates between the coefficients by using quad-trees. One of most efficient algorithms based on the EZW concept is the SPIHT algorithm. The SPIHT codec adds to the EZW algorithm a set partitioning technique. The SPIHT algorithm combines adaptive quantization of the wavelet coefficients with coding. The SPIHT algorithm presents a scheme for progressive coding of coefficient values when the most significant bits are coded first. This property allows control of the compression rate. The produced bitstream is further compressed losslessly by the application of adaptive arithmetic coding. In this invention, we can use either SPIHT or EZW.